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ABSTRACT 

We present a statistical framework which can be used to determine the contribution of 
an unresolved population of pulsars to the gamma-ray background. This formalism is 
based on the joint analysis of photon time series over extended regions of the sky. We 
demonstrate the robustness of this technique in controlled simulations of pulsar pop- 
ulations, and show that the Fermi Gamma-ray Space Telescope can be used to detect 
a pulsar contribution as small as 0.1% of the gamma-ray background. This technique 
is sensitive to pulsar populations with photon fluxes greater than ~ 1CP 10 cm~ 2 s _1 . 
The framework is extensible to arbitrarily complex searches for periodicity and can 
therefore be tailored to specific applications such as all-sky surveys and studies of the 
Galactic center and globular clusters. 
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1 INTRODUCTION 

The Large Area Telescope (LAT) onboard the Fermi 
Gamma-ray Space Telescope (FGST or Fermi) is a unique 
instrument that collects energetic photons from the whole 
sky, at an energy and spatial resolution as well as in an en- 
ergy range that offers a new window on high energy astro- 
physics. In almost three years since the launch of FGST, the 
sensitivity of the instrument has facilitated the discovery of 
new classes of objects, including gamma-ray pulsars. Over 80 
gamma-ray pulsars have been discovered in the Fermi— LAT 
all-sky data (see Abdo et al.| (|2010e|), and also Abdo et al. 
(2008J [ 2011] |2009b| |2010b| |2009c|); [Cognard et al.| (|20lTfI 



| Ray et ah (2011 1; Ransom et al. (2011 1; Saz Parkinson et al 



( 2010[ ); |Abdo et al.| ( |2010a|b[ ); |Keith et al.| ( |20li| ) 

While it is expected that more pulsars will be discov- 
ered as the baseline of the experiment is extended, most will 
remain undetected because their fluxes are below the sen- 
sitivity level of current detection techniques. These pulsars, 
as a population, contribute to the diffuse gamma-ray back- 
ground. Untangling the contributions to this background 
has been a subject of great interest, not only in the con- 
text of pulsar physics ( Watters fc Romani|2011 1 , but also in 
studies aimed at understanding the gamma-ray background 
near the Galactic center ( Buckley, Hooper fc Rosner||2011 



Hooper fc Goodenough|2011||Abazajian|2011||Dobler, Cho-| 
lis fc Weiner|2011||Boyarsky, Malyshev fc Ruchayskiy|2010[ 



Hooper fc Linden|2011 1, as well as a means to extract faint 
signals from exotic sources, such as dark matter (Malyshev, 



Cholis & Gclfand 2010) and antimatter (Gendelev, Profumo 
fc Dormody|2010 | 



In this paper we propose a new statistical search strat- 
egy that can be used to learn about the cumulative contri- 
bution of pulsars to the gamma-ray background. This tech- 
nique is an example of a general philosophy /strategy that we 
advocate, which is based on the concept that even though in- 
dividual data samples may not contain a detectable source, 
the statistics of a large number of samples contains informa- 
tion about the sources (see also |Geringer-Sameth fc Koushi-| 
|appas| ( f2010[ )). For the particular case we are studying here, 
even when a pulsar is not detected within a region of the 
sky, the data from that region will still contain informa- 
tion. When a large amount of such data is aggregated one 
can identify a statistical signature of the presence of pulsars 
even though the individual objects may not pass sensitivity 
thresholds. 

Such a statistical analysis can reveal the properties of 
the unresolved pulsar population. Application of this tech- 
nique to Fermi-LAT data can place bounds on the cumula- 
tive contribution of pulsars to the gamma-ray background 
which are independent of known sources. It is therefore a 
complimentary approach to the individual studies of bright 



pulsars with Fermi (|Buccheri, 
et al.|2006 l. 



Sacco & Ozel 1987 Atwood 



* email: alex_geringer-sameth@brown.edu 
f email: koushiappas@brown.edu 



We begin in Sec. [2] by describing the general strategy 
that can be used to learn about populations of objects when 
each individual one is undetectable on its own. We discuss 
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this in the context of the unresolved pulsar contribution to 
the gamma-ray background. In Sec. [3] we propose a spe- 
cific implementation involving the statistics of the maximum 
peaks in a collection of power spectra. It is developed in the 
framework of classical hypothesis testing, where the goal is 
to reject the null hypothesis that no pulsars are present in 
the gamma-ray sky. This includes the development of the 
statistical tests used to reject this null hypothesis. In Sec. [4] 
we make predictions for this method as applied to data from 
Fermi-LAT and show that under a wide range of circum- 
stances Fermi should be able to discover the presence of 
unresolved pulsars. Additionally, we show that individual, 
flux-unresolved, pulsars may be discovered based only on 
analysis of their time series. We discuss ways to extract the 
cumulative pulsar contribution to the background, which re- 
quires making assumptions about parameters of describing 
the pulsar population. Finally, in Sec. [5] we outline how this 
technique can be generalized to use more powerful tests for 
periodicity and discuss caveats which can affect the sensi- 
tivity of the method. 



2 GENERAL METHODOLOGY 

The detection of a pulsar at high significance relies on statis- 
tical tests performed on a collection of photon arrival times. 
(At radio frequencies, where the vast majority of pulsars 
have been discovered, the time series is in the radio inten- 
sity, not photon counts.) 

For the sake of simplicity assume that a certain statisti- 
cal test boils down the entire time series into a single num- 
ber, a "score'Q which is supposed to represent the "level 
of periodicity present" . The higher the number the stronger 
the periodic signal. "Detecting" a pulsar is an exercise in 
classical hypothesis testing and one needs to take into ac- 
count the fact that even if there is no pulsar present the 
score may be high because of random chance. Specifically, 
one needs the probability distribution for the score condi- 
tioned on the null hypothesis that there is no pulsar present. 
The question is asked, "What are the chances that the score 
would have been as high as measured if there was no un- 
derlying periodicity in the time series?" If the answer is, for 
example 0.3%, then a pulsar is said to be detected at 99.7% 
(or "3a") significance. In this example, the value of 0.3% is 
called the false alarm probability and in practice a 3a de- 
tection is hardly convincing. Usually, discoveries are claimed 
when the false alarm probability is less than 6 x 10 -7 , a "5a" 
detection threshold. 

The dominant factor in the detectability of a gamma- 
ray pulsar is the number of its photons which are collected 
by the LAT (i.e. the pulsar's photon flux). So far, Fer mi has 
detected pulsars with fluxes as low as 10 -8 cm -2 s" 1 ( Abdo 



et al. 2010el. These are pulsars whose time series are ex- 
tremely unlikely to have been generated by a non-periodic 
process — unlikely in the sense just discussed. However, it is 
quite likely that for every pulsar with such a flux the Galaxy 
contains a great many more with much smaller fluxes. If we 
assigned a periodicity score to the time series of these faint 

1 Throughout this article "score" is used in this sense and has 
nothing to do with the statistical concept of score defined as the 
derivative of the log-likelihood. 



pulsars the false alarm probabilities would be considerably 
greater. Most of them would be of order 1. Individually, 
these pulsations are undetectable with current data and pe- 
riodicity tests. 

However, what if one computes the periodicity score for 
40,000 time series, i.e. for every 1 square degree pixel on the 
sky? A few of these pixels will contain bright pulsars that 
will be unambiguously detected (these are the pulsars that 
are discovered using current pulsar search techniques). It is 
possible that many more pixels contain pulsars which are 
not obvious in the data (i.e. their periodicity scores are not 
improbably high), while most of the pixels will likely contain 
no pulsars at all. The goal then is to infer the presence of 
the undetected population of pulsars. 

The method we propose in this manuscript is based on 
a very simple observation: The periodicity scores from many 
separate time series, taken as collection, will be skewed to- 
ward larger values due to the presence of pulsars. By ana- 
lyzing the distribution of scores we can learn about a pop- 
ulation of objects whose individual members remain unde- 
tected. 

This general idea is not limited to the study of the 
galactic pulsar population. In fact, the concept of analyz- 
ing a collection of individually ambiguous signals to learn 
about a population underlies many studies of diffuse back- 
grounds. As an example, measuring the empirical counts 
PDF in sky pixels has been exploited in the study of blazars 



ter annihilation in substructure l 


Baxter & Dodelson 


2011 


Baxter et al.| 


2010| |Dodelson et al.||2009| |Lee, Ando &| 


Kamionkowski 


2009, Siegal-Gaskins|2008|), as well as pulsars 


(Faucher-Giguere & Loeb||2010 Siegal-Gaskins et al.|2010|. 



In these cases, the fact that the PDF differs from Poisson in- 
dicates that localized sources contribute to the background 
(even though any single "hot pixel" does not constitute a 
detection of an individual source.) 

A very simple example can illustrate the idea. Imagine 
we have a collection of 40,000 coins of which 98% are fair 
while the other 2% are rigged to land on heads 90% of the 
time. We get to flip each of the coins once and then try 
to answer the question, "Are there any unfair coins in this 
sample?" On the basis of one flip we have no way of say- 
ing whether any individual coin is fair or not. But perhaps 
the overall distribution of flip results can reveal information 
about the population of unfair coins. For example, suppose 
this experiment results in getting the expected number of 
heads: 40000 x (0.98 x 0.5 + 0.02 x 0.90) = 20320 heads. We 
pose the hypothesis test: if the coins were all fair what is the 
probability of getting 20320 or more heads? The answer is 



40000 /.fwwwA 

P (> 20320)= ]T 40 °°° (0.5) 40000 

i=20320 V / 



0.0007. (1) 



That is, there is a 0.07% chance of getting the results we did 
if every coin were fair. The hypothesis that all the coins are 
fair has been rejected with greater than 99.9% significance. 

Translating this scenario into pulsar language, each coin 
represents a one square degree patch of the sky. Flipping a 
coin corresponds to computing the periodicity score from 
that pixel's photon time series. Heads is a "high" score and 
tails a "low" one. If a pixel contains a pulsar the periodicity 
statistic gives a high score 90% of the time. The periodicity 
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score for a pixel with no pulsar present has equal chances 
of being high or low and one can not make any definitive 
claims based on the results of an individual measurement. 
However, the cumulative number of "high periodicity scores" 
from all 40,000 square degrees is strongly inconsistent with 
"no pulsars" . 



2.1 Cookbook 

The strategy discussed so far is general but can be decom- 
posed into several specific tasks. Here, we will outline the 
necessary steps, and in Sec. [4] we will develop a specific re- 
alization of this procedure which has been designed for ap- 
plication to Fermi-LAT data. 

The first step is to take the gamma-ray events in a re- 
gion of the sky and divide them into spatially separated time 
series. This can be done based on a simple pixelization of 
the sky or by collecting the photon time series from many 
promising locations (we will address these choices in Sec.[5|. 
Some preprocessing of the data should also be performed 
(e.g. applying a barycenter correction to each time series 
which corrects for the detector's motion with respect to the 
"fixed" solar system barycenter), as well as detector-specific 
corrections (e.g., see the Fermi Science Support CenteiQ. 

Next, a periodicity test statistic is chosen and applied 
to each time series. The choices for the test are numerous. 
We will detail a straightforward choice in Sec. [4] In general, 
the requirement is that one must assign a "score" to each 
time series which in some sense reflects the level of peri- 
odicity present. The test should be tailored to the type of 
objects one is searching for. For millisecond pulsars (MSPs), 
for example, it may not be necessary to take into account 
the effects of spin-down (see Sec. [5j|. 

It is essential to quantify the response of the test statis- 
tic to a white noise time series, i.e. an uncorrelated se- 
quence of photons which was not generated by a pulsar. 
Specifically, one needs the probability distribution for the 
score under the null hypothesis that no pulsar is present. 
This is called the null distribution. In the coin flipping 
example we used above, this probability distribution was 
Po (heads) = Po (tails) = 0.5. In some cases the null dis- 
tribution can be derived analytically. For more complicated 
periodicity tests the distribution can be found by simply 
running the periodicity test many times on randomly gen- 
erated white noise time series. 

Finally, given the collection of scores from the various 
time series, one tests the collection as a whole for deviation 
from the null distribution. There are a number of statistical 
tests that can be used for this purpose. Choices include the 
Kolmogorov-Smirnov and Anderson-Darling tests as well as 
the traditional \ 2 test of the binned histogram of scores. For 
the present application, we introduce an additional test, the 
A-test. It is designed to be sensitive to a very small tail of 
high periodicity scores (see next section and Appendix for 
more details). 



2 http: ll± ermi . gsf c .nasa. gov/ ssc/ 
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3 SPECIFIC IMPLEMENTATION 

In this section we present a methodology based on the above 
strategy. The goal is to detect the presence of unresolved 
pulsars by jointly examining the photon time series from 
numerous pixels in some area of the sky. For the sake of 
simplicity, we will assume that the pulsar period derivatives 
are very small. This particular implementation is appropri- 
ate for a search for the cumulative contribution of MSPs 
(Lorimer (20081, Ransom ( 2007 1), but can easily be gener- 



alized to the case where period derivatives are significant. 



3.1 Choice of periodicity test 

We need a numerical quantity, calculated from the measured 
photon data from each pixel on the sky, that describes the 
level of periodicity present in the time series. For this exer- 
cise the periodicity score of a time series is chosen to be the 
normalized peak magnitude of the power spectrum. We now 
explain what this quantity represents and how to compute 
it from a list of discrete photon arrival times. 

The Fourier transform is an alternate representation of 
the time series which highlights the various sinusoidal com- 
ponents that make up the signal. If a pulsar light curve is 
a pure sine wave its Fourier transform is a delta function 
spike at the pulse frequency. A well-used technique in pul- 
sar searches is to take the squared magnitude of the complex 
Fourier transform, called the power spectrum, and search 
for peaks in this function. The statistics of the power spec- 
trum for both random data (e.g 



Middleditch 



Ransom, Eikenberry & 



( 2002 i) and for data which contains a signal 



Groth (19751; Vaughan et al.| ( |1994[ ) have been well studied 
in general and in the context of pulsar searches. 

If photons arrive at times t\,ti, . . . , tjv we treat the sig- 
nal as a train of delta pulses at these times: 



3=1 

Plugging this into the definition of the continuous-time 
Fourier transform yields 



S(f) 



°? N 

/ e- 2mft s(t)dt = J2 e ~ Mft3 - 



(2) 



The unnormalized power spectrum is the absolute 
square magnitude of the Fourier transform. It is normalized 
by dividing by the mean power at each value of /. For data 
which contains systematic noise, calculating a running mean 
is required and may not be trivial. [Ransom, Eikenb erry fc| 
Middleditch ( 2002 ) present several techniques, including us- 
ing a running mean or a running median (divided by ln(2) 
) to normalize the power spectrum. For gamma-ray data at 
the high frequencies associated with MSPs there is likely 
no systematic non-white noise spectrum contaminating the 
time series. In this case (pure white noise) the mean is sim- 
ply equal to the number of discrete photon events in the 
time series. Therefore we search for peaks in the normalized 
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power spectrum P(f) defined as 
P(f) = ' 

cos(27ri/tj 




(3) 



+ 



^sin(27ri ftj) 
J=i 



We are only interested in the maximum of this quantity, and 
so computationally it is not necessary to store the entire 
Fourier transform in memory at any one time. This obvi- 
ates the need for the 10 billion point Fast Fourier Trans- 
forms (FFTs) that would be required for time series that 
are years long. Instead, one can calculate the power spec- 
trum by making incremental steps in the frequency, only 
saving the maximum power seen so far. This procedure is 
trivially parallelized by dividing the frequency interval to 
be searched into subintervals and searching each of these for 
its highest peak. |Ransom, Eikenberry fc Mid dlcditch ( |2002| ) 
provide trigonometric recurrences which can keep track of 
the the two sums in Eq. [4] as / is incremented in small steps 
without having to compute sines and cosines. 

The power spectrum is not an independent quantity 
for all values of /. It is a standard result from the study 
of discrete Fourier transforms that independent frequency 
"bins" have width 1/T, where T is the elapsed time over 
which the data was taken. For example, a three year LAT 
observation results in a width of 10 -8 Hz for each indepen- 
dent frequency bin. In searching for MSPs we would like to 
search over a frequency range corresponding to pulsar pe- 
riods between, say, 1 ms and 100 ms. In order to perform 
the search for peaks in the power spectrum we would first 
compute P(f) starting at / m j n = (100 ms) -1 = 10 Hz and 
then take steps of sizcr|<5/ = 1/T ~ 10 -8 Hz until reaching 



fn 



(1 ms) = 1000 Hz. Therefore, the exploration of 



the normalized power spectrum for each time series requires 
searching A/bins frequency bins, where 



A/bins = (/max - /min)T W 9 X 10* 



(4) 



In general, pulsar light curves are more complicated 
than sine waves which results in the Fourier transform hav- 
ing a series of spikes at integer multiples of the pulsar fre- 
quency. This fact motivates many pulsar searches to look 
for spikes in the sum of the first k harmonics of the power 
spectrum. Here we perform a more simple analysis that does 
not include the statistical details of searching the harmonic 
sum. However in practice, the pulsar search may be more 
sensitive if the highest harmonic-summed peak is used as 
the test statistic. We defer the discussion of various choices 
for the test statistic to a later section. 

In summary, we compute the normalized power spec- 
trum for the photon arrival time series for each pixel on the 
sky. The peak power in the power spectrum (in the frequency 
range of interest) is assigned to that pixel as its "periodic- 
ity score" . We will now explore the probability distributions 
describing the scores. 



3 In practice, one usually searches using a smaller step size in or- 
der to accurately explore each potential peak in the power spec- 
trum. However, this does not change the number of independent 
frequency bins searched. 



3.2 Statistics of the power spectrum peak for 
random data 

For each pixel the maximum of the power spectrum is a 
random variable. Following standard notation we call the 
random variable X. A specific realization (or measurement) 
of X is denoted by a lowercase x. If a pixel does not contain 
a pulsar, we assume that its power spectrum is just white 
noise, i.e. there are no periodic signals present in the fre- 
quency range of interest. In this case, the normalized power 
in each independent frequency bin is distributed according 
to an exponential distribution with a mean of 1 (e 



som, Eikenberry fc Middleditch (2002 1 ) 



Ran- 



Under the null hypothesis of no pulsars the score 
X is the maximum of A/bins independent exponentially 
distributed random variables. The cumulative distribution 
function (CDF) F(x) is the probability that all of the Abins 
random variables are less than x. This is simply equal to 
[Fx{x)] Mhi ™, where Fi(x) = 1 - exp(-z) is the CDF for a 
single exponentially distributed variable. The value of Abins 
is large (Eq. |4| and we can therefore make the following 
approximation, 

1 A/bins 



F(x) 



[1 



exp(- 



logAf binB ) 



1 - 



Ab in 



.<Vh 



(5) 

This result holds to high precision when Abins ~ 10 10 . 

The above expression shows that X is distributed ac- 
cording to what is known as a Gumbel distribution, some- 
times called an "extreme value distribution" . The probabil- 
ity distribution falls off extremely rapidly to the left of the 
mode at x = log A/bins and has a less steep tail to the right. 
Because log A/bins is a location parameter of the distribution 
the width of the Gumbel distribution does not change as 
A/bins increases. Also note that as the observation time in- 
creases the distribution shifts to the right at a logarithmic 
rate. This has important consequences that we discuss later. 
Looking ahead, as the observation time T increases, a pul- 
sar's power will grow in proportion to T while the random 
power it competes with grows only as logT. 

It is easy to invert F(x) to find 



x = logA"bins - log(- log-F). 



(6) 



Therefore, given a uniform deviate F between and 1, Eq.[6] 
can be used to transform it into a Gumbel distributed ran- 
dom variable. 



3.3 Statistics of the power spectrum peak when a 
pulsar is present 

The only distribution needed in order to perform an exper- 
iment that tests whether pulsars are present in the gamma- 
ray background is the null distribution given by Eq. [5] The 
test is simply whether the collection of time series is consis- 
tent with none of them containing any pulsar signal. In that 
case the score X for each time series is distributed as Eq. [5] 
However, in order to test the sensitivity of this method 
we need to be able to simulate situations where pulsars are 
present in the sky. In fact, to learn anything about the details 
of the pulsar population one needs some sort of model for the 
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way pulsars contribute to the background. Here we discuss 
how the presence of a pulsar affects the chosen periodicity 
statistic. We will return later to the question of extracting 
population parameters from the time series data. 

When a pulsar contributes photons to the time series, 
the peak of the power spectrum is distributed differently. 
In this case X is distributed as the maximum of two vari- 
ables. The first is a random variable representing the power 
in the bin at the pulsar's frequency. The second is a Gumbel 
distributed variable corresponding to the maximum power 
in the other (A/bins — 1) frequency bins. For frequency bins 
which are not at the pulsar's frequency, the pulsar photons 
contribute to the Fourier transform as if they were ran- 
domly distributed along with all the other photons. That 
is, the normalized power spectrum for the (A/bins — 1) other 
frequency bins is a white noise spectrum. We have already 
shown that the maximum power that will be found in these 
(A/bins — 1) bins is distributed according to F(x) (Eq.|5|. 

In order to determine the height of the normalized 
power spectrum for the bin at the pulsar's frequency we have 
to go back to the definition of the Fourier transforrrj^] The 
Fourier transform (Eq. [2| is seen to be the sum of unit vec- 
tors in the complex plane, one vector for each photon in the 
time series. In the case of white noise, each of these JV vectors 
has a random direction and the sum can be thought of as the 
endpoint of a random walk. This gives rise to the power in 
one frequency bin being distributed according to the expo- 
nential distribution with scale parameter JV. More precisely, 
let y be the sum of JV randomly directed 2-dimensional unit 
vectors. The direction of y will be uniformly distributed be- 
tween and 2-7T. The squared length of y will be distributed 
according to 



Prob« < jy| 2 < C + dC) 



-C/N 



N 



(7) 



It is easy to see that the normalized power in such a fre- 
quency bin, given by |y| 2 /JV, is exponentially distributed 
with scale parameter equal to 1, as stated above. 

Consider now a time series where JV S photons come from 
a pulsar and JVb g are uncorrelated background photons, such 
that the total number of photons is N — N s + JVb g . We ex- 
amine the Fourier bin at the pulsar's frequency and consider 
the idealized case where all the pulsar power lies in this sin- 
gle frequency bin with no power in harmonics. In this case 
each vector in the sum in Eq. [2] over the JV S pulsar pho- 
tons points in the same direction. It therefore has a length 
equal to JV S . The other JVb g background photons point in 
random directions and their sum in the Fourier transform is 
given by a randomly directed vector whose squared length 
I is distributed according to Eq. [7] with JV replaced by JVb g . 
To get the value of the normalized power spectrum for this 
frequency bin we take the squared length of the sum of the 
"signal vector" and the "background photon vector" and 
divide by the total number of photons in the time series. 
Defining P p to be the normalized power in the frequency 
bin at the pulsar's frequency we have 

Pp = jV i Ns2 + 1 + 2iVs ^ cos(f9) ] ' 



4 This paragraph is based on the geometric interpretation given 
in Vaughan et al. 1 19941. 



The power spectrum height is seen to be a random variable: 
the quantity I is distributed as I ~ (l/JVb g ) exp(— Z/JVb g ) and 
9 is a uniform random variable between and 2-7T. 
We introduce the following new variables: 



fb 



JV 3 . 
VJV ' 

. Nbz 
" JV 



yjvTTJV, 



JV s + JV b 



(8) 



(9) 



The first can be thought of as a signal to noise term repre- 
senting how many photons in a pixel are due to a pulsar vs. 
background. The second measures the fraction of photons in 
a pixel which are not due to the pulsar. In terms of these 
the normalized power spectrum becomes 

P p = aS 2 + I' + VT> cos 0, (10) 

where the variables 6 and I' are distributed according to 

I' ~ i e - ! '//b 
fb 

6 ~ Uniform [0, 2tt]. 

In Eq. |10| a is introduced to take into account the frac- 
tion of the pulsar's total power that lies in this single fre- 
quency bin. If the light curve of the pulsar were a perfect 
sine wave all of the signal power would lie in the bin at the 
fundamental frequency and a — 1. In more realistic situ- 
ations the power will be divided up into higher harmonics 
and a may be less than 1. 

For reference, we note that the probability distribution 
of P p has been worked out analytically in Groth (19751, 



which also contains general results that may be of use when 
considering more complicated tests for periodicity. In partic- 
ular, the probability distribution for the sum of an arbitrary 
number of harmonics in the power spectrum is also derived. 



3.4 Rejecting the null hypothesis of "No pulsars" 

As described above, each sky pixel is assigned a periodicity 
score X which is defined to be the peak height of its normal- 
ized power spectrum. The goal is to take this collection of 
X values and perform a statistical test of the following null 
hypothesis: The time series for every pixel is nothing but 
white noise, i.e. no pulsars are present in any of the pixels. 
More precisely, we ask if the collection of measured X values 
is consistent with each score being drawn from the distribu- 
tion in Eq. [5] (i.e. generated by a random white noise time 
series) . 

This task can be accomplished by a number of statistical 
methods. Here, we use a new test developed specifically for 
this application. In this section we outline how the test works 
and refer the reader to the Appendix for details. 

It is desirable (and possible) to use a classical hypothesis 
test to learn about the sensitivity of this method. The idea 
is to boil the collection of measured X values into a single 
test statistic we call A. The quantity A should, in some 
sense, indicate the overall level of periodicity present in the 
gamma-ray sky, just as X did for a single pixel. Small values 
of A should indicate "less evidence for periodicity" than do 
large values of A. 
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The " A test" is based on the quantity (see Appendix), 



A = 



^-log[l-F(xi)] 



- NlogN + logM 



(11) 

where the Xi are the measured scores (normalized power 
spectrum peaks) for each of the N time series and F(x) is 
the CDF of the null distribution given by Eq. [5] The test is 
designed to give more weight to time series with large scores. 

The test statistic is treated as a random variable and 
its probability distribution under the null hypothesis (that 
every sky pixel contained only non-periodic, random pho- 
tons) is quantified. A significance threshold is chosen and 
the critical value A* is defined so that if the null hypothesis 
holds, then the probability that A < A* equals the cho- 
sen significance. For example, if we want to perform a "3er" 
search one finds A* such that P(A < A*) = 0.997. If we 
find that the observed value of A is in fact greater than A* 
the null hypothesis is to be rejected at "3cr" significance. In 
other words, it would be extremely unlikely to measure such 
a high value of A if there were no pulsars. This indicates 
that pulsars contribute to the gamma-ray background. 



4 APPLICATION TO FERMI-LAT 

We now turn to the question of detecting the presence of 
pulsars in the gamma-ray sky using current data. We assess 
the conditions where the proposed formalism is successful in 
rejecting the null hypothesis of "no pulsars" in the diffuse 
background as measured by the LAT instrument on board 
Fermi. In this section, we will demonstrate the robustness 
of this method by generating simulations which contain a 
controlled population of pulsars with known properties. We 
utilize the maximum normalized power periodicity test along 
with the A test as described above. 

Assume that a region of the sky is isotropically popu- 
lated with pulsars that all have the same flux, $ p , defined 
as photons per area per time in some energy range. These 
pulsars contribute a fraction 7 of all the photons received 
by the LAT in this energy range. That is, of all the photons 
that LAT detects over the entire sky a fraction 7 of these 
originated from pulsars each having a flux 4>p. The projected 
number density of pulsars is given by a p (number of pulsars 
per solid angle). 

The average flux the LAT measures is given by F to t in 
units of photons per area per time per solid angle (in the 
relevant energy range). In addition to pulsars we assume a 
uniform, isotropic background flux Fbg (same units as Ftot). 
The independent parameters of this model are $ p and 7. 
The background flux is chosen to make up the difference 
between the pulsar contribution and the observed total flux. 
Specifically, 



and 



Ftot = Fb e + o-p$p, 



7-FtOt = 0><I>j; 



(12) 



(13) 



These two equations determine Fb g and a p in terms of & p , 
7, and the observed Ftot. These equations are more easily 
interpreted by multiplying through by the solid angle of the 
survey and by the observation time and effective area of the 



detector. Then F tot becomes the total number of photons re- 
ceived by the LAT over the entire survey area, a p becomes 
the total number of pulsars in the survey area, & p the num- 
ber of photons received from each pulsar, and F Dg the total 
number of background (non-pulsar) photons received over 
the survey area. Solving the above equations we find 



and 



F hs = (1 - 7)i?U, 



7^t< 



(14) 



(15) 



We assume a value of F to 



8.72x10" 



in the energy range [0.8-6.4]GeV ( |Abdo et al.|2010d| |. This 
includes the energy range in which pulsars are most impor- 
tant relative to the total flux ( Abdo et al.|2010e |. 

In order to generate simulated data, we need a survey 
area and pixel size. We choose the pixel size, 0, to be I 
square degree, and we will use two choices for the survey 
area: 40,000 square degrees which represents the all-sky sur- 
vey, and 1,000 square degrees, which roughly represents the 
inner ~ 32 x 32 degrees around a region such as the Galactic 
center. 

We must evaluate Eqs.[8]and[9]to generate a normalized 
power for each pixel that contains a pulsar. The number of 
background photons in a pixel is 



jV bg = F bg A eff T, 



(16) 



where A c a is the (orbit-averaged) effective area of LAT (2000 
cm 2 ) and T is the observation time (3 years). The number 
of pulsar (signal) photons in a pixel which contains a pulsar 
is 

As = $p AeS T. (17) 
Inserting these quantities in Eqs.[8]&[9j we have 



F bg + $ p 



$>vA cS T 



-1/2 



(18) 



7 + 



F tt 



F tot 



VF tot 0A ctt T, 



and 



A = 



l+(%/F tot 0)- 



(19) 



For a given choice of $ p and 7 we can use these last 
two equations along with Eq. [10] to generate a normalized 
power in a pixel that contains a pulsar [^] For simplicity the 
simulations were performed using a = 1. Consequences of 
relaxing this assumption will be discussed later. 

We explore the parameter space to see when pulsars 



5 There are many choices for <E> p and 7 that give a number of 
pulsars which is larger than the number of pixels, i.e. a p > 1. 
When this is the case we need to generate a normalized power 
for each pulsar in the pixel, a peak power from the other ~ A/"bm S 
frequency bins and then take the maximum of all these to be the 
periodicity score X for the pixel. We have found that for the range 
of parameter space we discuss the extra pulsars in each pixel do 
not change the results. Therefore, we run the simulations with 
at most one pulsar per pixel (though cr p is allowed to be greater 
than 1). 
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Figure 1. A demonstration of the statistical power of the method 
to detect the presence of pulsars over the entire sky. The color cod- 
ing represents the probability of rejecting the null hypothesis of 
"no pulsars" at 99.7% significance. <3? p is the photon flux of each 
individual pulsar in the energy range [0.8 — 6.4] GeV. The quan- 
tity 7 represents the fraction of the total gamma-ray background 
due to pulsars. Solid contours give the number density of pulsars 
(in units of pulsars per square degree). The proposed method can 
reveal the presence of a pulsar population contributing as little 
1CT 3 of the diffuse gamma-ray background. Note that, within 
the range of pulsar fluxes shown, every individual pulsar is /Jus- 
unresolved because <J> P is less than LAT's point source sensitivity 
threshold. Many of these flux-unresolved sources may be individ- 
ually discovered based solely on an analysis of their time series: 
the dashed line represents the 5cr detection threshold for individ- 
ual pulsars based on the height of their power spectrum peak (see 
text for details). 



will be detected by this method using the A statistic defined 
above. We choose a value of A* corresponding to a 99.7% 
( "3cr" ) detection. For each pair of values $ p and 7 we create 
1,000 realizations. For each realization we simulate 40,000 
(all-sky) and 1,000 (Galactic center) values of X (one for 
each pixel) and compute the A statistic. Out of the 1,000 
trials we count the number in which the null hypothesis is 
rejected. The fraction of trials in which the null hypothesis 
is rejected is the sensitivity (or power) of the proposed test. 
For example, if for a particular choice of $ p and 7 we find 
that in 900 out of 1,000 simulations the null hypothesis is 
rejected (i.e. A > A* in 900 of the simulations), then there 
is an 90% chance of making a "3cr" detection of the presence 
of pulsars. 



4.1 Results 

Figure [l] shows the results of the parameter space scan over 
values of $ p below 10 -9 cm 2 s _1 and over the full range of 
7 from 10 -5 to 1 for a simulated all-sky survey of 40,000 
square degrees. The colour-coding corresponds to the power 
of this method to reject the null hypothesis that there are no 
pulsars at 99.7% ("3<r") significance. In the dark red region 
the null hypothesis is practically guaranteed to be rejected. 
In the blue region the null hypothesis will be rejected only 



0.3% of the time (as expected for a 99.7% significance thresh- 
old). The solid contours correspond to the number density of 
pulsars (in units of pulsars per square degree) as computed 
using Eq. [15] 

There are two competing factors which shape the tran- 
sition between the sensitive and insensitive regions of pa- 
rameter space. The plateau at small values of 7 represents 
the limit of low numbers of pulsars. Obviously, if there are no 
pulsars in the sky there is no signal to be detected. Within 
the flux range explored here the A test is not sensitive if 
there are fewer than ~ 10 pulsars in the 40,000 pixels. 

The vertical transition is explained by the fact that pul- 
sars must contribute the highest peak in the power spectrum 
in order to be detected by the periodicity test. As the flux of 
each pulsar is increased (moving to the right in Fig. [T]) the 
power spectrum peak at the pulsar's frequency will eventu- 
ally become the highest peak in the power spectrum. This 
then causes the non-Gumbel-ness of the pixel scores which 
is detected by the A test. 

We can view this as a requirement that the quantity 

? A/bins 



Pp (Eq. 101 be comparable to log A/bins, the mode of the 
distribution for the maximum normal ized power in the case 
of no pulsars. The aS 2 term in Eq. 10 is most important 



term in Eq. 

in governing the transition. Because the Gumbel distribtion 
only contains a location parameter we can write an approx- 
imate equation describing the vertical part of the sensitivity 
transition: 



aS 2 ~ logA/" bins . 



(20) 



The left hand side is an estimate of the height of the peak 
corresponding to the actual pulsar signal. The right hand 
side is the maximum power in the other A/bins — 1 frequency 
bins. Only when the left hand side is greater than the right 
hand side will the method be able to reject the null hypoth- 
esis of no pulsars. This is because the periodicity statistic 
we have chosen is not sensitive to pulsar peaks which are 
subdominant in the power spectrum. 

The photon fluxes of individual pulsars in the simulated 
parameter space are all below the point source sensitivity of 
the LAT ( jAtwood et al.|2009[ ). The pulsars in the simulation 
would be undetected by Fermi as bright sources. Therefore, 
"blind searches" would not consider these pulsars as candi- 
dates for periodicity searches. Such objects truly contribute 
to the diffuse background. 

Nevertheless, if we measure the power spectrum for a 
pixel which contains an unresolved pulsar the spike at the 
pulsar's frequency might be large enough to constitute a de- 
tection of a periodic source. To estimate when this occurs we 
consider the following hypothesis test. We measure a power 
spectrum peak height of x, and ask "What is the probabil- 
ity that at least one peak in any of the observed time series 
has a value of x or higher if there were no periodic sources 
present in the data?" The answer again follows from the 
Gumbel distribution (Eq. [5]l but with A/bins replaced with 
A/bins x A/" p ix, i.e. the number of independent frequency bins 
for each time series multiplied by A/pix, the number of time 
series considered (in this case 40,000). The quantity F(x) is 
the significance of this peak. 

In the region to the right of the dashed line in Fig. [l] 
individual pulsars would be detected at 5a based on the 
height of the power spectrum peak derived from their pixel's 
time series. The region's shape is governed by an equation 
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Figure 2. Same as Fig. [l] but for an observation area of 1,000 
square degrees corresponding to a study of the Galactic center. 
Here, solid contours depict the total number of pulsars present 
in the observed region. The dashed line denotes the 5<r detection 
threshold of individual pulsars based on power spectrum peak 
height as in Fig. [I] 



similar to Eq. |20| except that the right-hand side is replaced 
by a peak height xs a such that 1 — F{x^ a ) — 5.7 x 1CP 7 , 
corresponding to a 5a detection. This suggests that simply 
computing the power spectra for the entire sky may turn up 
detections of pulsars which are too faint to be flux-resolved. 

In Fig. [2] we show the results of a similar simulation 
but for an observation of only 1,000 square degrees of sky. 
This situation represents a study of the galactic center, a 



region whose source population is of great interest (Buck 
ley, Hooper fc Rosner|[20TT| [Hooper fc Goodenough poTTf 
Abazajian 2011 Dobler, Cholis fc Weiner||2011| Boyarsky, 



Malyshev fc Ruchayskiy|2010||Hooper fc Linden|2011[ ). The 
shape of the sensitivity region is similar to the all-sky sur- 
vey. The stripe pattern is caused by the discrete addition of 
pulsars to the survey area. The method is sensitive to the 
presence of pulsars even when individual pulsars are unre- 
solved in both flux and Fourier power. The sensitive (dark 
red) region is larger in the all-sky survey than in the galactic 
center study. This demonstrates that the statistical test ben- 
efits from larger numbers of measured time series (assuming 
equal fluxes and number densities of pulsars). 



5 DISCUSSION 

Using the approximation to the sensitivity transition given 
by Eq. [20]we can estimate how the result of the simulations 
discussed in the previous section will scale with changing pa- 
rameters. The Gumbel distribution we have been exploring 
has the beneficial property that the location parameter goes 
as the logarithm of the number of independent trials. We 
expect this to be a general feature of any periodicity test. 
Thus, as observation time T increases the right hand side of 
Eq. |20| increases as logT. The left-hand side increases in pro- 
portion to & p 2 A c ffT/0. Thus, this technique benefits from 
longer observation times, larger effective areas, and smaller 



pixel sizes (i.e. future gamma-ray observatories) in the same 
"root N" way that conventional searches do. 

The main difficulty in the outlined strategy lies in choos- 
ing a good test of periodicity and in the computational chal- 
lenge of computing it many times for the different time se- 
ries. Traditionally, pulsars are searched for either by taking 
a Fourier transform of the time series or by folding the time 
series in the time domain at many different trial periods. 

Regular pulsars do not have constant periodicity but ex- 
perience spin-down (magnetic braking) and glitches. These 
complicating factors force statistical tests for periodicity to 
be performed on a large grid of trial frequency derivatives 
or on short stretches of the time series ( |Atwood et al.|2006| 
Buccheri, Sacco fc Ozel|1987 1 . In Fourier space the changing 
period of the pulsar acts to spread its signal power over many 
frequency bins, diluting the peak amplitude. Millisecond 
pulsars, on the other hand, have extremely stable r otations, 
with period derivatives on the order of 10~ 19 s/s (Lorimer 



|2008[ ). Even over observation periods of years the frequency 
of many MSPs will not drift into neighboring Fourier bins 
( Ransom|2007 1. Thus, the number of trials performed when 
computing the test statistic can be significantly lower than 
for regular pulsars. 

Additionally, MSPs are thought to form in binary star 
systems. Because binary systems are more common than sin- 
gle stars most galactic pulsars are likely members of a binary 
pair that have been spun up into MSPs ( Alpar et al.||1982| 



Radhakris hnan fc Srinivasan|1 982 Bha ttacharya fc van den| 



Heuvel 1991 Tauris & van den Heuvel 2006 1 In addition, 



it has recently been suggested that MSPs might dominate 
normal pulsars in their contribution to the gamma-ray back- 
ground ( Faucher-Giguere fc Loeb|2010 1 . Millisecond pulsars 
are also older, have had more orbital trips around the galaxy, 
and therefore are more likely to be found at higher galactic 
latitudes than normal pulsars. Therefore, these pulsars may 
be important contributors to the so-called "extra-galactic" 



or isotropic gamma-ray background (Siegal-Gaskins et al 



20101 



5.1 Caveats and Improvements 

We have been optimistic in some areas and overly simplistic 
in others. Here we review some of the practical difficulties 
in performing this test on LAT data and point out the sim- 
plifications we have made and how they affect the results. 

The most obvious difficulty we have glossed over is the 
fact that pulsar light curves are not simple sine waves. This 
has the effect of dispersing the power in the frequency bin 
centered at the pulsar's frequency into higher harmonics. 
The simple periodicity test we proposed (maximum normal- 
ized power) is almost certainly not optimal for the case when 
power is found at higher harmonic frequencies (see below for 
ways to try to recover this power). We have left room in the 
analysis (see Eq. 10 1 for a reduction of a, designed to ac- 
count for the effect of power being dispersed into other fre- 
quency bins. Equation [20] suggests that the pulsar flux one 
is sensitive to goes as (since S scales proportionally 

to $ p ). 

While millisecond pulsars are extremely stable and do 
not experience glitches or suffer from rapid spin-down, their 
rotation is not completely constant. It is therefore probable 
that some power is dispersed into neighboring frequency bins 
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by non-negligible period derivatives. The techniques used to 
try to recover this power involve performing many analyses 
with different trial period derivatives. Specifically, the arrival 
times of the photons are corrected to account for a spin- 
down effect and then the periodicity search is performed on 
this modified time series. The decrease in sensitivity due to 
spin-down increases as the observation time increases. 

As millisecond pulsars are found in binary systems, the 
orbital motion of the binary can cause distortions in the 
Fourier spectrum of the time series. Essentially, the orbit 
of the pulsar causes a doppler shift in its period which dis- 
perses Fourier power into different frequency bins. Methods 
have been proposed that can sweep up this power ( |Ransoni^ 
Cordes fe Eikenberryp003 1 . Such methods can be incorpo- 
rated into a more advanced periodicity test. 

Errors in source position are known to affect the de- 
tectability of individual pulsars. The first step in analyzing 
a time series is to correct for motion of the detector with re- 
spect to the pulsar. This correction depends on an accurate 
"barycentering" procedure, which in turn relies on precise 
knowledge of the pulsar's position. In searching the back- 
ground for unresolved pulsars we have no information as to 
where the pulsars are located within the pixels, and this 
affects the quality of the barycentering. 

The most important consideration that goes into a real- 
istic application of the proposed method is the choice of peri- 
odicity statistic. In practice, one is bound by finite computa- 
tional resources — ideally, one would perform a detailed time 
series analysis on every pixel in the sky, including searching 
over trial periods, period derivatives, and other ephemera. 
We have been simplistic in the choice of the maximum nor- 
malized power as a test statistic. A first generalization is 
to search harmonic sums of the normalized power spectrum. 
This would take into account that pulsar light curves are not 
sine waves. Considering the harmonic sum of the power spec- 
trum is an attempt to recover the as much signal power as 
possible. The statistics of such a test are relatively straight- 
forward to compute. 

In addition, there are several choices of tests for peri- 
odicity currently in use to search for pulsars in radio data 



and in gamma-rays . The H test (de Jager, Raubenheimer 
fc Swanepo"el||l989[ ) and the Z\ test (Buccheri et al.||1983| ) 



are based on binning the photon arrival times by phase for 
a given trial pulsar period and then checking whether the 
distribution of phases is consistent with random. These tests 
require a guess for the pulsar period. However, it is computa- 
tionally intensive to calculate the statistic for every possible 
value of the pulsar period for a large sample of pixels. More 
recently, a time-differencing technique ( Atwood et al.|2006 1 
has been proposed to overcome some of these computational 
challenges and has been very successful in discovering new 
pulsars with Fermi-LAT ( |Abdo et al.|2009a| |20T0e] ). 

To adapt these tests to the present task, we propose 
to first find the power spectrum of the time series (or its 
harmonic sum) and take the n highest peaks to be trial 
periods for the more advanced algorithms. The number of 
trials n would need to be adjusted based on computational 
resources and the choice could be calibrated by examining 
the power spectra from time series which are known to con- 
tain gamma-ray pulsars. An automated analysis pipeline can 
be conceived in which one would perform a cursory scan of 
the time series looking for semi-significant peaks and then 



perform additional, computationally intensive scans of these 
peaks, assigning a periodicity score at the end. 

Besides computational cost one has to balance two fac- 
tors when deciding on a periodicity test. The test should be 
as sensitive as possible to the presence of periodic signals 
but should also minimize the number of "trials". A large 
number of trials raises the possibility that a random signal, 
by chance, could appear periodic. In our case the number of 
trials was the number of independent Fourier bins that were 
scanned when looking for peaks. As the number of "trials" 
grows it is more likely to find a random outlier that mimics 
periodicity. 

Any periodicity test or analysis procedure can be 
adapted to the search for unresolved pulsars. The key in- 
gredient is the null distribution of the periodicity scores. 
For example, an arbitrarily complex analysis pipeline can 
be established which takes a time series and outputs a peri- 
odicity score. The inner- workings of the pipeline can involve 
scanning over trial periods and period derivatives. It can in- 
clude identifying promising peaks for more careful scanning. 
Once the procedure is set, one simply runs it many times 
on uncorrelated photon time series (i.e. white noise). The 
resulting set of periodicity scores constitutes the null dis- 
tribution. The pipeline is then applied to actual measured 
time series and the resulting scores are collectively checked 
for inconsistency with the null distribution. 

We can illustrate the effect of different tests using the 
sensitivity plot in Fig. [T] Different periodicity tests would 
shift both the sensitive (red) region and dashed line together. 
However, the scaling is not necessarily fixed. The dark red 
region of parameter space to the left of the dashed line re- 
mains the most interesting. It is here where searches for 
individual pulsars would fail but where the collective statis- 
tics would succeed in revealing their presence. The size and 
shape of this region likely depends on which periodicity test 
is chosen. We plan to explore other tests in future work. 

Furthermore, the division of a region of the sky into 
spatially separated time series (step one in Section 2.1) can 
also be optimized. Instead of breaking the sky into pixels 
and taking the time series of each one, an alternate tech- 
nique is to only search promising sky locations for evidence 
of periodicity. One could consider only "bright spots" or 
"hot pixels", regions of the sky with a signal to noise ratio 
greater than 1, say. Alternatively, the candidate locations 
can be chosen from lists of known sources (see Fermi bright 



source list, Abdo et al. (2009dl, Abdo et al. (2010c l), or from 



pulsar candidate locations in blind searches. The later have 
been previously analyzed for pulsations but have not been 
jointly searched for unresolved pulsars. These strategies have 
several advantages. The computational burden would be re- 
duced because of the fewer number of time series to scan. 
The barycenter correction would be improved by the better 
localization of the sources' positions. A priori, hot pixels have 
a higher chance of containing pulsars than randomly selected 
pixels, leading to a larger fraction of the searched pixels that 
contain pulsars (effectively increasing a p in Fig. [TJ . 

Because the analysis is sensitive only to the highest 
power spectrum peak it is almost completely insensitive to 
the possibility that there may be multiple pulsars contribut- 
ing to a single time series. However, this situation likely oc- 
curs in globular clusters and in the galactic center region, 
both places conceivably containing important populations 
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of pulsars. A periodicity statistic should be tailored specif- 
ically to studies of these regions. A simple generalization 
of the periodicity test would be to take the top n high- 
est peaks in each time series instead of just the highest. 
Then we would have n periodicity scores from each pixel 
instead of one. Alternatively, one could count the number 
of peaks with height greater than some threshold. The score 
from each pixel would be this integer number. (In both cases 
the search could take place using the harmonically summed 
power spectra.) 



5.2 Pulsar population parameter estimation 

This analysis begs the followup question of how we can learn 
the details of the pulsar population from studies like this, 
where individual pulsars remain undiscovered. In particu- 
lar, it is of great interest to determine what fraction of the 
gamma-ray background is due to unresolved pulsars (the 
value of the quantity 7 in the simulations of Sec. |4j|. The 
detailed extraction of population parameters from the col- 
lection of periodicity scores requires some kind of model- 
ing of the population. However, we can use the simplified 
model presented here to place interesting constraints on the 
number of pulsars with certain fluxes without any detailed 
modeling. 

In the simulations of Sec. [4] we assumed that every pul- 
sar had the same flux. This is obviously false if we claim 
that the simulated pulsars make up all the pulsars in the 
sky. However, the simulated pulsars can instead be inter- 
preted as a "slice" of the number function of pulsars. 

An important description of the pulsar population is the 
number density of pulsars 0" P ($ P ) with flux greater than <I> P . 
This function can be used to define the total contribution 
from pulsars: 



1 

-ftot 



(21) 



The simple simulations of Sec. [4] can be used to con- 
strain (jp($ p ) as follows. Imagine that we have performed 
a test over the whole sky (Fig. [T]) but failed to reject the 
null hypothesis. At each flux <I> P we can draw a line straight 
upwards in Fig. [l] until we reach the transition to the dark 
red region. Let the number density of pulsars simulated at 
this transition point be given by a p (& p ). Then we can claim 
that the true number density function at this flux <J P (& P ) 
must be less than cr p (<l? p ). If this were not the case then 
there is a 99.7% (in this example) chance that the statisti- 
cal test would have detected the presence of these pulsars. 
This constraint relies on the choice of a and in practice the 
choice should be calibrated using known pulsar light curves. 

If we are willing to make some assumptions about the 
shape of ct p ($ p ) and only allow its overall normalization to 
vary we can make stronger statements. In this case we could 
actually simulate a population of pulsars for different choices 
of normalization and find the sensitivity of the method to 
each choice. The test will become sensitive above some crit- 
ical value of the normalization. Depending on whether the 
test rejects or does not reject the null hypothesis we could 
then place a lower or upper bound on the normalization of 
the number density function. This bound would then imme- 
diately translate into a bound on the total contribution of 



pulsars to the background (Eq. |21[ ). There are several mo- 
tivated choices for the shape of cr p ($ p ) which depend on 
the spatial distribution of pulsars jFaucher-Gi guere fc Loeb] 
(2010). In reality, however, the population of gamma-ray 
pulsars is completely u nconstrained at flu xes below about 
10~ 8 photons cm" 2 s _1 jAbdo et al.|2010e |. 

In addition, one can analyze the measured distribution 
of periodicity scores using conventional \ 2 minimization. In 
this case it is necessary to know what the distribution of 
scores will be as a function of the pulsar population param- 
eters. One then can bin the measured scores and find the 
best fitting population parameters. The pulsar population 
models can be made as complicated as one likes — the anal- 
ysis requires a scan over this parameter space looking for 
regions whose score distribution matches the observed one. 
We defer applications of these techniques to the LAT data 
in future work. 



6 CONCLUSIONS 

In this manuscript we propose a new technique whose appli- 
cation to Fermi-LAT data can reveal the extent to which pul- 
sars contribute to the gamma-ray background. The method 
is based on the cumulative statistics of photon time series 
that are binned spatially. The motivation behind this ap- 
proach lies in the general idea that even though individual 
pulsar searches may be unsuccessful, information from unde- 
tected pulsars is still measurably encoded in the gamma-ray 
background. 

In general, current pulsar searches are based on the ev- 
idence of a source at a particular location. These sources 
are subjected to a battery of periodicity tests, and careful 
analysis of LAT data has already revealed the presence of 
gamma-ray pulsars. However, it is likely that large num- 
bers of pulsars are beyond the current reach of LAT to even 
identify their associated events. These pulsars (with very 
weak signals) will contribute to the diffuse gamma-ray back- 
ground. 

Our main results are: 

• The proposed technique has the ability to discover a 
pulsar contribution to the gamma-ray background if the 
fraction due to pulsars is greater than 10 -3 . 

• It is sensitive to a population of pulsars whose individ- 
ual photon fluxes are as low as 10~ 10 cm -2 s _1 . 

• Using the photon time series derived from a specific 
location on the sky, one can discover individual pulsars with 
photon fluxes down to about 6 x 10 -10 cm -2 s _1 , which is 
below the current point source sensitivity threshold. 

• By considering only "hot pixels" or current blind 
search candidates the sensitivity of the method is increased 
markedly. 

• Any periodicity test or analysis pipeline can be applied 
to the search for the unresolved population. The only re- 
quirement is the response of the test to uncorrelated photon 
time series. This allows the technique to be optimized for 
any given application (e.g. all-sky surveys, galactic center, 
globular clusters, etc.). 



The method proposed in this work takes advantage of 
all events in the diffuse gamma-ray background and gives in- 
formation about the population of unresolved pulsars. The 
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importance of this task goes beyond pulsar astrophysics. It is 
manifestly apparent that a detailed understanding of astro- 
physical backgrounds is vital in any gamma-ray observation, 
including surveys of astrophysical sources (e.g., blazars), as 
well as studies of more exotic and hypothetical contribu- 
tions (e.g., annihilating dark matter). It is therefore of ex- 
treme interest to apply this technique to current and future 
gamma-ray data. 

We acknowledge useful conversations with Jacqueline 
Chen, Ian Dell'Antonio, Andrew Favaloro, Scott Field, 
Richard Gaitskell, Elizabeth Hayes, Kyle Helson, Julie 
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Jason Su. SMK thanks the Aspen Center of Physics for hos- 
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APPENDIX A: THE A TEST 

In this appendix we provide details about the statistical test 
we used in this paper. The test is designed to determine if 
a collection of observations is inconsistent with having been 
drawn from a given null distribution. It is meant to be sen- 
sitive to a small upper tail in excess of what is predicted by 
the null distribution. Although motivated by the application 
to pulsars the A test has nothing to do with astrophysics and 
may be used in any statistical study. 



Al Motivation 

Recall the situation presented in the text. We have a col- 
lection of periodicity scores (denoted Xi) and want to test 
whether the collection is consistent with having been drawn 
from the null distribution (in this case a Gumbel distribu- 
tion) . The goal is to boil the collection of scores down into a 
single number A and then study the distribution of A under 
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the null hypothesis. The quantity A is meant to reflect the 
overall level of periodicity in the sample. 

The critical value A* is defined by the property that, 
if the null hypothesis is true, the probability that A is less 
than A* is e.g., 99.7%. To be precise, A is a function of 
the collection Xi. If the x's are each drawn from the null 
distribution then the probability that A is less than A* is 
0.997, or whatever the desired significance is. 

Different choices of A may be more or less powerful. In 
general, the power of a test is defined as the probability that 
the null hypothesis is rejected when the null hypothesis is, 
in fact, false. If it is unlikely that A is above some critical 
value A* even when there are many pulsars present in the 
sky a poor definition for A has been chosen. Unfortunately, 
only in special, simple cases is there a "uniformly most pow- 
erful" test. In the particular case we are studying here there 
are many degrees of freedom associated with the alterna- 
tive hypothesis. For example, the light curves of pulsars and 
their number density are functions which must be specified 
by many parameters. As a result there is no uniformly most 



powerful test in this case. (See e.g. Stuart & Ord (19871 for 
an more detailed discussion.) 

In order to choose a powerful statistical test we must ex- 
amine the behavior of the collection of x's in the case where 
pulsars are present. Consider a pixel which contains a pul- 
sar. The only way the a;-value of this pixel will contain any 
information about the pulsar is if the peak in the normalized 
power spectrum is actually due to the pulsar. Under the null 
hypothesis, each x is drawn from the Gumbel distribution 
in Eq. [5] The effect of pulsars is to skew the distribution 
towards higher values of x: the pixels with a pulsar have a 
chance of replacing the peak power in a random frequency 
bin with the power at the pulsar's frequency. Based on these 
considerations we would like to choose a statistical test that 
puts more weight on higher x values. 

There are a wide variety of statistical tests that are 
in common use. The Kolmogorov-Smirnov (KS) statistic 
is commonly used in astronomy. Kuiper's extension of the 
KS statistic gives more weight to the tails of the distribu- 
tion. This would be beneficial for looking for the excess in 
large a;-values. The Anderson-Darling statistic is used more 
rarely but also gives extra weight to the tails. Likelihood ra- 
tio statistics are another option, though these require some 
knowledge of the alternative hypothesis that one is testing 
for. It is known that likelihood ratio tests are the most pow- 



erful tests for "point" hypotheses Stuart & Ord ( 1987 1. They 



are based on the likelihood function for the data under var- 
ious hypotheses, and should therefore exploit all the infor- 
mation available in the data. 

The proposed A test statistic is designed to be sensitive 
to the upper tail of a distribution. It shares properties with 
the Anderson-Darling and KS tests and can also be inter- 
preted as a likelihood-ratio test. Unlike these other tests, 
however, the distribution of the A test statistic under the 
null hypothesis is very simple (a gamma distribution). It is 
expected to be powerful (like a likelihood test) but also very 
easy to use (no sorting of the data and no lookup tables) . 



if this collection is consistent with being drawn from a given 
probability distribution (the null distribution). Below, this 
collection of numbers will also be referred to as the "data" 
or the "samples" . 

When looking for an extended tail in a collection of mea- 
sured quantities we noticed that it is often useful to look at 
the logarithm of the empirical survival function (SF) of the 
data. The empirical SF is defined as 1 — Fn(x), where Fn(x) 
is the empirical cumulative distribution function (CDF). 
Simply put, the SF at some value x is the fraction of the 
sample values which are greater than x. Thus, at x = — oo 
the empirical SF equals 1 and decreases by 1/N every time 
x crosses one of the measured values, where TV is the sample 
size. This empirical SF can be compared to the theoretical 
SF for the case where the data come from the null distribu- 
tion. For the null distribution, the survival function is simply 
1 — F(x), where F(x) is the usual cumulative distribution 
function for the null distribution. 

When comparing the logarithm of the empirical and 
theoretical SFs any excess at large values of x becomes more 
pronounced, even if only a small fraction of the samples 
are at such large values. Therefore, we order the data by 
increasing a;- value and define the A statistic as 

1 N 

A =-/^J2 { lo S I 1 - - log [1 - F( Xi )]} , (Al) 

* i — 1 

where x\ < X2 < ■ ■ ■ < xn and Frf(xi) — (i — l)/N, the 
empirical CDF. We can make some simplifications to the 
first term in the sum as : 



^og[l -F N ( Xi )] = ^log 



(i-l) 



N 



log 



TV- 1 
TV 



= log [TV (TV — 1) (TV — 2) • • • 1] — log [TV^] 
= log TV! -AT log AT. (A2) 
Inserting this back into the definition of A we have 



A = 



- TV log TV + log TV! 



(A3) 

The statistics of A is governed by the term in curly brack- 
ets. In this sum the numerical ordering of the rr's does not 
matter since the sum is over all of them. The distribution 
of A under the null hypothesis is now straightforward to 
find. For any random variable X with CDF F the quantity 
F(X) is distributed uniformly in the interval between and 
1. This implies that 1 — F(X) is also uniformly distributed 
on this interval. Now, the negative logarithm of such a uni- 
formly distributed variable is distributed according to the 
exponential distribution with scale factor 1. Therefore, un- 
der the null hypothesis the quantity 



G=£)-]og[l- J(a*)] 



(A4) 



A2 Details 

In this subsection we present the details of the A test. The 
task is to take a collection of numerical values and determine 



is the sum of N exponentially distributed random variates. 
This sum is described by the well-known gamma distribution 
(also called the Erlang distribution in this case) with shape 
parameter N. The inverse CDF of the gamma function then 



© 0000 RAS, MNRAS 000, 000-000 



Pulsars in the diffuse gamma-ray background 13 



provides the critical value A*. For instance, to find the value 
of A* under which there is a 99.7% chance of measuring A 
(under the null hypothesis) one determines the value of G* 
that satisfies 



0.997 = 



G * x N ~ 



(N-iy. 



dx. 



(A5) 



The quantity G* is then inserted into Eq. |A3| replacing the 
term in curly brackets. The resulting value of A is A*. If 



for a given sample of N x-values the quantity A (Eq. A3 1 is 



greater than A* then one can reject that the sample came 
from the distribution with CDF F(x) at 99.7% significance. 

A3 Properties of A 

Of course, there is no reason to include the constant terms in 



Eq. A3 One can just take the test statistic to be G (Eq. A4 1, 
the only quantity that depends on the data. Then G* , dis- 
cussed above, is the critical value for the test statistic. (In 
fact, this is how we actually performed the simulations.) 



However, the definition we have given for A (Eq. A3 1 has 



a nice asymptotic property for large sample sizes (i.e. as 
N — > oo). The central limit theorem says that the gamma 
distribution converges to a normal distribution with mean 
N and standard deviation V^V- In the same limit the con- 
stant term Eq. |A2| converges to — N as can be seen using 
the approximation for log(AM) found in every statistical me- 
chanics textbook (e.g. Reif (1965]), section A. 6). Therefore 
as N — > oo the distribution for A converges to a standard 
normal distribution (i.e. normal with mean and variance 
!)• 

The A test statistic is similar to the KS and Anderson- 
Darling statistics in that is based on the CDF of the null 
distribution. The CDF has the nice property that it is dis- 
tributed uniformly (if the null hypothesis is true). This al- 
lows the null distributions for the KS, Anderson-Darling, 
and A test statistics to be found analytically. 

The specific application of the A test statistic shown in 
this paper can also be interpreted as a likelihood-ratio test. 
The null distribution is given by the Gumbel distribution 
with a peak at log A/bins- Imagine that the alternative dis- 
tribution for the x's follows the null distribution for values 
of x less than log A/bins but does not fall off for higher val- 
ues. This is supposed to represent the situation when pulsars 
are present: there are more large values of x. The likelihood 
ratio is the ratio of the alternative PDF to the null PDF 
(as functions of x). When this quantity is large it indicates 
that the alternative describes the sample better than the 
null does. The likelihood ratio is the product of these ratios 
for each Xi. It is usually easier to work with the logarithm 
of this quantity which is the sum of the logarithms of the 
individual likelihood ratio terms. 

Let us see how each term in the log-likelihood ratio 
compares to each term in the G statistic (i.e. each term 
in the curly bracketed sum in Eq. |A3[ ). If x is less than 
logA/bins both statistics contribute approximately 0. In the 
case of the likelihood ratio this is because the null and al- 
ternative PDFs are defined to be the same there (so the 
log of their ratio is 0). It is also easy to see from Eq.|5]that 
when x is less than log Abms, F(x) is close to 0. If x is greater 
than log A/bins the quantity 1 — F(x) becomes approximately 
exp(-(a;-logAbins)) and so -log(l-F(a;)) ~ a; — log A/bins- 



For the likelihood ratio when x > log Abms the alternative 
hypothesis PDF is 1 and the null PDF is approximately 
exp(— (x — log A/bins)). Thus the logarithm of this ratio is 
also approximately x — log Abms- 

For all values of x, therefore, the A statistic (based on 
the quantity G) behaves just like a likelihood ratio test that 
is designed to pick up an extended upper tail in the sample. 
This implies that the A test should be a powerful test in 
looking for such a tail. Moreover, the null distribution of A 
has a particularly simple form (a shifted and scaled gamma 
distribution) and converges to the standard normal distri- 
bution when the sample size is large, making A an attractive 
addition to the current library of tests. 
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